% This program is written using Dyanare v.5.4, which needs to be downloaded from www.dynare.org. To run the program, this mod file as well as the steady state file Tables2and3_steadystate should be in the same folder. 
%The steady state is derived from the solution to the equilibrium along the lines described in Online Appendix C4. Further computation details arre given in Section 3.   
%The structure of the program follows the standard Dynare format. The notation of the variables and parameters follows that in the paper. The program generates impulse response functions (IRFs), a subset of which are utilized in Figures 1 and 2.
%Following usual practice, variable names are sometimes altered slightly to avoid confusion with Matlab functions. Thus, alph=alpha, gamm=gamma, sigmaes=sigma. Q is an auxiliary variable, and x1, x2 and sigmaf are defined for computational convenience. 

%close all;
%Endogenous variables

var p1 p2 U Q EOmega I1 I2  K1 K2 A1 A2 D1 D2 W P u1 I1tilde S1tilde S1 S2 c1 c2 C  IR2 Rf rf  R1 R2 Er1 Er2 Erp1 Erp2 Y1 Y2 gy1 gy2 Rf_real rf_real  c IR1 ir1 ir2 gc w gw gir1 gir2 pi gpi;
predetermined_variables K1 K2;

%Exogenous variables
varexo e1 e2;  

%Parameters
parameters alph, rho1, rho2, gamm, h1, h2, m1, m2, delta1, delta2, phi1, phi2, x1,x2,sigmaf, psik1,psik2, varphi1,varphi2,sigmaes,N1,N2,thet, eta;

%Calibration;
h1=0.05;
h2=0.05;
m1=0.1;
m2=0.1;
delta1 = 0.05;
delta2=0.05;
phi1 = 0.5;
phi2=0.5;
gamm = 5;
eta=0.525;
psik1=0.49;
psik2=0.49;
alph=0.99;
sigmaes=4.8;
varphi1=8;
varphi2=8;
N1=7;
N2=54;
rho1=0.935;
rho2=0.935;
sigmaf=(sigmaes-1)/sigmaes;
x1=(sigmaes*(1-psik1)+psik1)/sigmaes;
x2=(sigmaes*(1-psik2)+psik2)/sigmaes;
thet=(1-gamm)/(1-eta);

model;
p1=W^(1/sigmaes)*P^sigmaf*phi1*(N1*A1*(u1*K1)^psik1)^(-1/sigmaes);
p2=W^(1/sigmaes)*P^sigmaf*phi2*(N2*A2*K2^psik2)^(-1/sigmaes);
D1=(p1-h1)*(A1*(u1*K1)^psik1)-(I1+0.5*varphi1*(I1/K1-delta1)^2*K1)-m1*K1;
D2=(p2-h2)*(A2*K2^psik2)-(I2+0.5*varphi2*(I2/K2-delta2)^2*K2)-m2*K2;
c1=(W/P)*(P*phi1/p1)^sigmaes;
c2=(W/P)*(P*phi2/p2)^sigmaes;
C=(phi1*c1^((sigmaes-1)/sigmaes)+phi2*c2^((sigmaes-1)/sigmaes))^(sigmaes/(sigmaes-1));
U= ((1-alph)*C^((1-gamm)/thet) + alph*Q^(1/thet))^(thet/(1-gamm));
Q=U(+1)^(1-gamm);
P=(phi1^sigmaes*p1^(1-sigmaes)+phi2^sigmaes*p2^(1-sigmaes))^(1/(1-sigmaes));
EOmega=(P/P(+1))*alph*(C(+1)/C)^((1-gamm)/thet)*C/C(+1)*U(+1)^(eta-gamm)*Q^((1/thet)-1);
W=N1*D1+N2*D2;
I1tilde=(-0.29*K1+12.61*K2);
S1tilde=(0.03-9.16*K1+48408.76*K2); 
(1+varphi1*(I1/K1-delta1))*(p1*(1-u1^psik1)+u1^psik1*sigmaes*(p1-h1))/p1=EOmega*((W(+1)^(1/sigmaes)*P(+1)^sigmaf*phi1*psik1*sigmaf*A1(+1)^sigmaf*u1(+1)^(psik1*sigmaf)*K1(+1)^(-x1)*N1^(-1/sigmaes)-psik1*h1*A1(+1)*u1(+1)^psik1*K1(+1)^(psik1-1))*(1-u1^psik1/p1*(p1-sigmaes*(p1-h1))*u1(+1)^(-psik1*sigmaf))-psik1*u1^psik1/p1*(p1-sigmaes*(p1-h1))*u1(+1)^(psik1/sigmaes)*(1-u1(+1)^(-psik1/sigmaes))*h1*A1(+1)*K1(+1)^(psik1-1)+(p1*(1-u1^psik1)+u1^psik1*sigmaes*(p1-h1))/p1*((1-delta1)*(1-varphi1*delta1)-m1)+0.5*varphi1*(((I1(+1)/K1(+1))^2-u1^psik1/p1*(I1tilde(+1)/K1(+1))^2-delta1^2*(p1*(1-u1^psik1)+u1^psik1*sigmaes*(p1-h1))/p1))+(1-delta1)/K1(+1)*varphi1*(I1(+1)*(p1(+1)*(1-u1(+1)^psik1)+u1(+1)^psik1*sigmaes*(p1(+1)-h1))/p1(+1)+I1tilde(+1)*((p1(+1)*(1-u1(+1)^psik1)+u1(+1)^psik1*sigmaes*(p1(+1)-h1))/p1(+1)-(p1*(1-u1^psik1)+u1^psik1*sigmaes*(p1-h1))/p1)));
(1+varphi2*(I2/K2-delta2))=EOmega*(psik2*W(+1)^(1/sigmaes)*P(+1)^sigmaf*phi2*A2(+1)^sigmaf*K2(+1)^(-x2)*N2^(-(sigmaes+1)/sigmaes)*(N2*sigmaes-1)/sigmaes-psik2*h2*A2(+1)*K2(+1)^(psik2-1)-m2+0.5*varphi2*((I2(+1)/K2(+1))^2-delta2^2)+(1-delta2)*(1+varphi2*(I2(+1)/K2(+1)-delta2)));
(p1-h1)*A1*K1^psik1*(1-u1^psik1) + S1tilde =S1;
S1=EOmega*(D1(+1)+S1(+1));
S2=EOmega*(D2(+1)+S2(+1));
R1=((W(+1)^(1/sigmaes)*P(+1)^sigmaf*phi1*psik1*sigmaf*A1(+1)^sigmaf*u1(+1)^(psik1*sigmaf)*K1(+1)^(-x1)*N1^(-1/sigmaes)-psik1*h1*A1(+1)*u1(+1)^psik1*K1(+1)^(psik1-1))*(1-u1^psik1/p1*(p1-sigmaes*(p1-h1))*u1(+1)^(-psik1*sigmaf))-psik1*u1^psik1/p1*(p1-sigmaes*(p1-h1))*u1(+1)^(psik1/sigmaes)*(1-u1(+1)^(-psik1/sigmaes))*h1*A1(+1)*K1(+1)^(psik1-1)+(p1*(1-u1^psik1)+u1^psik1*sigmaes*(p1-h1))/p1*((1-delta1)*(1-varphi1*delta1)-m1)+0.5*varphi1*(((I1(+1)/K1(+1))^2-u1^psik1/p1*(I1tilde(+1)/K1(+1))^2-delta1^2*(p1*(1-u1^psik1)+u1^psik1*sigmaes*(p1-h1))/p1))+(1-delta1)/K1(+1)*varphi1*(I1(+1)*(p1(+1)*(1-u1(+1)^psik1)+u1(+1)^psik1*sigmaes*(p1(+1)-h1))/p1(+1)+I1tilde(+1)*((p1(+1)*(1-u1(+1)^psik1)+u1(+1)^psik1*sigmaes*(p1(+1)-h1))/p1(+1)-(p1*(1-u1^psik1)+u1^psik1*sigmaes*(p1-h1))/p1)))/((1+varphi1*(I1/K1-delta1))*(p1*(1-u1^psik1)+u1^psik1*sigmaes*(p1-h1))/p1);
R2=(psik2*W(+1)^(1/sigmaes)*P(+1)^sigmaf*phi2*A2(+1)^sigmaf*K2(+1)^(-x2)*N2^(-(sigmaes+1)/sigmaes)*(N2*sigmaes-1)/sigmaes-psik2*h2*A2(+1)*K2(+1)^(psik2-1)-m2+0.5*varphi2*((I2(+1)/K2(+1))^2-delta2^2)+(1-delta2)*(1+varphi2*(I2(+1)/K2(+1)-delta2)))/(1+varphi2*(I2/K2-delta2));
Rf=(1/EOmega);
Rf_real=Rf*(P/P(+1));
rf=Rf-1;
rf_real=Rf_real-1;
Er1=R1-1;
Er2=R2-1;
Erp1=Er1-rf;
Erp2=Er2-rf;
c=log(C);
gc=c-c(-1);
IR1=I1/K1;
IR2=I2/K2;
ir1=log(IR1);
gir1=ir1-ir1(-1);
ir2=log(IR2);
gir2=ir2-ir2(-1);
Y1=p1*A1*(u1*K1)^psik1;
Y2=p2*A2*K2^psik2;
gy1=log(Y1)-log(Y1(-1));
gy2=log(Y2)-log(Y2(-1));
w=log(W);
gw=w-w(-1);
pi=log(P);
gpi=pi-pi(-1);
K1(+1)=K1*(1-delta1)+I1;
K2(+1)=K2*(1-delta2)+I2;
log(A1)=rho1*log(A1(-1))+e1;
log(A2)=rho2*log(A2(-1))+e2;
end;
steady(nocheck);
shocks;
var e1;stderr 0.03;
var e2;stderr 0.03;
corr e1, e2 = 0.0;
end;
stoch_simul(order=2, pruning, periods=590000, irf=10) S1 S1tilde S2 u1 p1 p2 D1 IR1 IR2 D2 W P C Er1 Er2 Erp1 Erp2 gc gw gpi gir1 gir2 Y1 Y2 gy1 gy2 rf_real; 




